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Abstract. The recovery type error estimators introduced by Zienkiewicz and Zhu use a 
recovered stress field evaluated from the Finite Element (FE) solution. Their accuracy 
depends on the quality of the recovered field. In this sense, accurate results are obtained 

00 

using recovery procedures based on the Superconvergent Patch Recovery technique (SPR). 

These error estimators can be easily implemented and provide accurate estimates. An- 

other important feature is that the recovered solution is of a better quality than the FE 

• • solution and can therefore be used as an enhanced solution. 

> 
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We have developed an SPR-type recovery technique that considers equilibrium and dis- 
cs 

placements constraints to obtain a very accurate recovered displacements field from which 
a recovered stress field can also be evaluated. We propose the use of these recovered fields 
as the standard output of the FE code instead of the raw FE solution. Techniques to 
quantify the error of the recovered solution are therefore needed. 

In this report we present an error estimation technique that accurately evaluates the 
error of the recovered solution both at global and local levels in the FEM and XFEM 



1 



1 Introduction 



frameworks. We have also developed an /i-adaptive mesh refinement strategy based on 
the error of the recovered solution. As the converge rate of the error of the recovered 
solution is higher than that of the FE one, the computational cost required to obtain a 
solution with a prescribed accuracy is smaller than for traditional /i-adaptive processes. 

KEY WORDS: error estimation; recovery; SPR; XFEM; error control; mesh adaptivity 



1 Introduction 



Here we briefly discuss the possibility of estimating the error in the recovered solution cr* 
evaluated from the FE solution cr h as part of the process of error estimation in energy 
norm based on recovery-techniques. Recovery error estimation techniques are based on 
the assumption that the recovered solution is more accurate than the FE solution. A 
sufficiently accurate estimation of the error in energy norm of the recovered solution could 
lead to refinement processes based on the accuracy of the recovered solution instead of 
that of the FE solution. This could result in a considerable reduction of the computational 
cost that would be particularly interesting, for example, in optimization processes whose 
efficiency would be significantly increased as the required level of accuracy would be gained 
with a considerable lower number of degrees of freedom. 

In [1] and the Ph.D. thesis [ ] we developed techniques whose objective was that of 
obtaining upper bounds of the error in energy norm by means of the use of recovery-based 
error estimators. This research yielded an expression of an error bounding technique 
which includes two types of terms. The first term is the expression of the ZZ error 
estimator. The second term is a correction term used to account for the lack of equilibrium 
of the recovered stress field (both in the interior of the domain and along the domain 
boundaries). Numerical experiments have shown that slightly manipulating the expression 
of the correction term one can estimate the error of the recovered solution, not only at a 
global level but also at element level. 

In this report we will show the initial developments that finally led to an empirical ex- 
pression based on this correction term that can be used to efficiently estimate the error 
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2 Model problem 



in energy norm of the recovered solution u*, initially evaluated to estimate the error of 
the FE solution u h . So, the report is structured following our developments along time. 
Section 2 shows the model problem. Section 3 shows the technique used to evaluate up- 
per bounds of the error based on nearly equilibrated recovered solutions that leads to 
the correction term previously mentioned. Section 4 describes the technique to obtain 
the recovered displacement and stress fields (u* and cr*) that will be used to estimate the 
discretization error in energy norm of the finite element solution u h . Section 5 will present 
an expression that could be used to evaluate the error in energy norm of the recovered 
solution and the empirical approximation to this expression that can be used in practice 
to efficiently estimate this error. Finally, the numerical examples showing the potential 
of the technique proposed for the estimation of the error in the recovered solution are 
presented in Section 6. 



2 Model problem 



Let us consider the 2D linear elastic problem. The unknown displacement field u, taking 
values in f2 C 1R 2 , is the solution of the boundary value problem given by 

-V<r(u) = b infi (1) 

a (u) • n = t onr w (2) 

u = onr fl (3) 



where and denote the Neumann and Dirichlet boundaries with dQ = Tn U Td and 
Tjy fl T D = 0. The Dirichlet boundary condition in (3) is taken homogeneous for the sake 
of simplicity. 

The weak form of the problem reads: Find u e V such that 

a(u,v)=Z(v) VvGV, (4) 
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3 Upper error bound of the FE solution from nearly equilibrated recovered solutions 



where V is the standard test space for the elasticity problem such that V = {v | v G 
H\{1),v\ Td (x) = 0}, and 

a(u,v):= f cr(u) : e(v)dfi = / cr(u) : S : cr(v)c/fi (5) 
Jn Jn 

/(v) := I b • vrffi + / t • vdT, (6) 

where cr and £ denote the stresses and strains, and S is the compliance tensor. 

The bilinear form a(-, •) can also be expressed in terms of stresses by formally introducing 
a(-, •) such that a(<r,r) := L<r: D -1 r<if2. Note that a(u, v) = a(<T,r). 

Let u h be a finite element approximation to u. The solution lies in a functional space 
V h C V associated with a mesh of finite elements of characteristic size h, and it is such 
that 

a(u h ,v) = /(v) VveV h (7) 

We focus on assessing the error e = u — u h . The error is measured in the energy norm, 
induced by a(-, ■). The quantity to be assessed is ||e|| 2 = o(e, e). 

3 Upper error bound of the FE solution from nearly 
equilibrated recovered solutions 

This section will explain with detail the basis of the proposed upper error bounding 
method. The ideas were presented in [1] and the Ph.D. thesis [2]. 
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3 Upper error bound of the FE solution from nearly equilibrated recovered solutions 



Proposition. Let o"* be a recovered field such that cr* n is continuous almost everywhere 
along any interior curve rcU (being n the outward normal vector to T) and 

—V • cr* = b + s a.e. in Q (8a) 

tr* • n = t + r a.e. on T N (8b) 

Note that if s = and r = 0, cr* is statically admissible. For small s and r (when 
compared to b and t, respectively) cr* is said to be nearly statically admissible, up to the 
equilibrium default terms. 

Then, the following expression holds: 

a(ff*, ct(v)) = Z(v) + [vsdn+[ v-rrfr Vv G V, (9) 

Jn Jv N 

and, as a direct consequence assuming that the Dirichlet boundary conditions are homo- 
geneous 

||u|| 2 = a(cr(u),cr(u)) < a(u*,u*) -2 / u-sdtt-2 / u-rdT (10) 

Jn Jr N 



Moreover, the error approximation c* = c* — cr(u h ) could also lead to an upper bound 
of the error norm: 



\\e\\ 2 = a(e, e) = a(cr e , cr e ) < a«, a*) - 2 / e • s dQ - 2 / e • r dT (11) 

Jn Jr N 

and this latter property stands also for non-homogeneous Dirichlet boundary conditions, 
up to oscillation terms. 

Proof. The expression in (9) is derived from (8) using the standard weighted residuals 
technique combined with integration by parts (and the fact that the test function v 



UPV 



5 



3 Upper error bound of the FE solution from nearly equilibrated recovered solutions 



vanishes on Tn) 

/ v • (b + s) dtt = - [ v • V • cr* dtl 
Jn Jn 

= - [ V • (cr* ■ v) dtt + [ Vv : cr* dtt 
Jn Jn 

= - I n- (cr* • v)dr+ I -(Vv + V T v) • a* dQ 
Jdn Jn 2 

= - / v • (cr* ■ n) dT + a(cr*, cr(v)) 
Jr N 

= - [ v- (t + r)rfr + a(cr*,tr(v)) 
Jr N 

Then (9) follows using the definition of Z(v). Thus, the proof of (10) is straightforward 
considering in (9) v = u (this requires that in the original problem, the Dirichlet boundary 
conditions states that u = on T D ) 

a(a*, <t(u)) = Z(u) + / u-sdtt+ / u-rrfr 

Jn Jr N 

= a(a(u),a(u))+ u-sdtt+ u-rdT (12) 
Jn Jr N 

and therefore 

< a(a(u) - cr*, <r(u) - cr*) =a(a(u), cr(u)) - 2a(cr*, cr(u)) + a(cr*, cr*) 

= — a(cr(u), cr(u)) — 2 / u ■ s dtt 

Jn 

-2 I u-rrfr + a(cr*,cr*). 
Jv N 



The proof of (11) is similar, by considering in (9) v = e (this can be done, up to os- 
cillation terms, even if in the original problem the Dirichlet boundary conditions are 
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4 Displacements and stress recovery 



non- homogeneous). Thus, 

a((T*, a(e)) = 1(e) + / e-sc(n+ / e • r dT = a(cr(e), cr(e)) 

+a(cr(w // ), <r(e)) + e-sdtt+ e • r dT 
Jq Jr N 

that results in 

a(cr*, cr(e)) = a(cr(e), tr(e)) + e-sd£l+ e-rdT. 

Jq Jr N 

Note that the following error representation has been used 

a(cr(e), cr(e)) = /(e) — a(a(u H ), tr(e)). 



Thus, (11) is proved using the same idea as before by simply considering the positiveness 
of a(cr(e) - v* e ,v{e) - a*). 



4 Displacements and stress recovery 



The evaluation of the correction terms in (11) involves the evaluation of the equilibrium 
defaults s and t and the exact error in the displacements field e. The evaluation of s 
and t depends on the recovery technique used to obtain cr* and we are able to derive 
expressions for the exact evaluations of these terms. However, e = u — u h can only be 
estimated. We estimate e as e m u* — u h . 

We have derived two alternatives for the evaluation of e. The first one consists in using 
a sequence of N meshes for the problem. The solution of the last mesh of the sequence 
would be considered as an enhanced solution u* such that in mesh i 

e (i) = u - uf i} w u h (N) - ufo, % = 1, N (13) 



UPV 



7 



4 Displacements and stress recovery 



/ e ■ sdQ < 2 


/ e • sdf2 


Jn 


Jn 



In [1] we evaluated the correction term for the lack of internal equilibrium considering 

<2|e| L Js| L2 (14) 



As previously mentioned, we are able to evaluate s® and therefore |s(i)|^. Then we 
use (13) to evaluate |e(i)| L for i — 1, ...,N whereas the value of | e (7V)| L2 is evaluated by 
extrapolation from the value of |e(jv-i)| L ■ 

This procedure was enhanced in [3] by evaluating the correction terms ES® for i = 1, N 

as 

ES® = -2 I e w ■ s (i) dfi « -2 / {u\ N) - u h (l) ) ■ s (i) dfi, z = 1, N - 1 (15) 
Jn Jn 

and then ES(n) by extrapolation from ES^n-i)- 

The two versions of this first alternative have a high computational cost as i) it requires 
the use of a sequence of meshes and ii)it requires projecting the information from mesh z 
to mesh N. Note also that the results obtained for the first meshes are quite accurate but 
this accuracy decreases in the last meshes, especially in the last mesh where the results 
have been obtained using extrapolation techniques. 



The second alternative consists in using a recovery technique that would directly provide 
an enhanced displacements field u* that could then be used to estimate the error in the 
displacements field as e u* — u h without requiring a sequence of meshes. In this section 
we will describe the key ideas of the recovery technique used to obtain the error estimates. 
The technique is also based on the SPR scheme [ I] , but considering the displacements field 
and the local satisfaction of the imposed displacement, boundary tractions and internal 
equilibrium. 

We define each of the components of the recovered displacements on the support (patch) 
of a vertex node % as a polynomial surface expressed as u* = paj where p is a polynomial 
basis and are the unknown coefficients. The degree of p is one order higher than the 
degree of the shape functions used for the interpolation of displacements. We use the 
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Langrange Multipliers technique to force u* to satisfy the Dirichlet boundary conditions 
if the patch contains the Dirichlet boundary. Furthermore, using Lagrange Multipliers 
over the coefficients a it is possible to enforce restrictions on <J*(x, y) = DLu*(x,y) so 
that u* satisfies the following equations: 

• Internal equilibrium equation. 

• Boundary equilibrium equation. 

• Compatibility equation (by default). 

The displacements in each element are evaluated from each of the vertex nodes of the 
element. To enforce continuity of the recovered displacements we use a partition of unity 
approach, defined by Belytschko as the Conjoint Polynomial enhancement [5]: 



Where N? represent the shape function of node i corresponding to the linear version of 
the element (note that patches are only formed around vertex nodes) and n v refers to the 
number of vertex nodes. Equation (16) provides a kinematically admissible displacement 
field u* evaluated as a function of local contributions u* from each of the different patches. 
We will use subindex u to indicate a kinematically admissible solution whereas subindex 
a will be used to indicate a solution with a continuous stress field. 




(16) 



i=i 
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With (17), we can evaluate the stresses cr* at each point as 







dx 

o I 

ay 

JL A. 
-dy dx- 



<j:(x,y)=-DLj2NUx,y)u*(x,y) 

i=i 

Tlu 

V) = DLN ^ vK(s, y) + N?(x, y)DLu*(x, y) 



(17) 



i=l 
s 



i=l 



discontinuous 



0%(x,y) continuous 



The expression to evaluate cr* has two parts, one discontinuous and another discontinu- 
ous. Taking only the continuous part we will have a continuous stress field cr* which is 
nearly-statically admissible as it is built using local contributions where the satisfaction 
of equilibrium has been considered. 

We refer to cr* as nearly-statically admissible field because due to the enforcement of 
continuity, a lack of internal equilibrium s is introduced as shown in (18) 



<(x,y) = ]TN^,y)DLu*(z,y) 

n v 

V<(:r, y) = V N *(^ VpLu^x, y) 
i=i 

Tlu Ylu 

V<(s, y) = Yl VN *^> y)DLu*(x, y) + N?(x, y) VDLu*(x, y) 



(18) 



i=l 



i=l 



V(T*=-b 



V<(x,y) = -s-b 
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5 Error norm representation for the recovered solution 



Moreover, there could also be a lack of contour equilibrium r over r^v, r = cr* • n — t, 
where n is the outward normal vector, if the polynomial used to describe the recovered 
field is not able to represent the applied boundary tractions t and/or if the boundary is 
not a straight line segment. 

5 Error norm representation for the recovered solu- 
tion 

At the same time we worked on an expression which allow us to evaluate the error of 
the nearly-statically recovered stress field cr* . Recall that this is a continuous stress field 
but not fully equilibrated, where s and r represent the lack of internal and boundary 
equilibrium, respectively. 

The recovery procedure yields a postprocessed stress field cr* which is taken as an en- 
hanced approximation to the exact stresses, cr(u), much more accurate than a(u h ). The 
recovered stress cr* has an equilibrium default 

s = -V-<-b (19) 

Pedro Diez proved that the following identity holds to evaluate the error in the recovered 
field: 

E* = — f s ■ e*jn = || 2 = ||u - <|| 2 = a ((cr - <), (cr - <)) (20) 
Jn 

being e* the error in displacements corresponding to the recovered solution cr* such that 
cr* = cr(u*) = cr(u — e*). Note that both u* and e* are not explicitly evaluated. 

Proof. 



UPV 



11 
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IdQ 



[ s-e* a dn= [ (V-< + b)-e* ff c 
Jn Jn 

= [ (V • cr«) - V • cr(u)) • e *dSl 
Jn 

= [ (-V • cr(e* a )) • e* a dn 
Jn 

= [ cr«) : e(e* a )dn - / < • a(e*J • ndT 
Jn Jan 

= [ cr«) : e(e* a )dn 
Jn 

= \K\\ 2 (21) 



Note that it has been assumed that 



/ e* a • <r(e* a ) ■ ndT = (22) 
Jan 



If we also consider this term we would have: 



lie* || 2 = - / s • e*jn + [ el- a«) • ndT (23) 
Jn Jan 

Considering cr(e*) ■ n = cr(u — u*) ■ n = cr(u) ■ n — cr(u*) • n = t — cr* • n and then 
considering (8b) we have cr(e*) = — r. Therefore ||e* || 2 can be evaluated as 



E* = \\e*J 2 = - I s • eldtt - I r • e*JT (24) 

an 



Operating with (20) considering C-S inequality we will have 



£f = - / e* a sdn<K\ L2 \s\ L2 (25) 
Jn 



UPV 



12 



5 Error norm representation for the recovered solution 



As u* would be a recovered displacements field, laying in a so-called 'broken space' richer 
than V h , we could assume that the L 2 norm of the recovered solution is smaller than the 
L 2 norm of the FE solution: 



le*l < lei 



(26) 



We would then obtain the following practical expression to evaluate an upper bound of 

1 1 2 

||e* || taking into account that the error of the FE solution can be estimated as e m e es 



with e P 



u; - u 



< |e* I r Isl r < lei T |s| T 



~ u;-u" 



\L 2 > S >L2 



\ e es\L 2 \ S \l 2 



E 



UB 



(27) 



Note that in (27) we have finally replaced |e*| La by |e es | L2 , i.e. we have replaced the error 
in the recovered solution u* by the estimated error of the FE solution u h , to obtain a 
bound of the error in the recovered solution. Following this idea we replace e* = u — u* 
by e es = u* — u h in (24), and defined the error indicator E* in (28a) to check if it could 
provide an indication of the error level in energy norm of the recovered solution cr* . We 
also defined the error indicators E\ and E\ as described in (28b) and (28c) to force the 
result to be positive: 



E\ = 



/ s ■ e es dQ — l 
Jq Jt 



Nel 



*5 = E 



e: 



r ■ e P „dr 



+ 



r • e es dT 



rnfii 



/ |s • e es | dVt + / |r-e es |dr 
Jn Jv 



(28a) 
(28b) 
(28c) 



where Nel is the number of elements of the FE discretization and Qj represents the domain 
of element i. 
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In (28b) the value of the integrals at each element are forced to be positive. In (28c) the 
integrands themselves are forced to be positive. Note that this is a reasonable assumption. 
For example, if we assume r = (see (21)) 

< a(e* ff ) £ (e*J = -s • e* a (=|s-e* ff |) (29) 

As s and e* are consistent (s would be the defaults of equilibrium corresponding to u*, 
whose associated error is e*) then < — s ■ e* . However, in (28c) e* has been substituted 
by e es . The terms s and e es are non-consistent and as a result — s • e es could be negative. 
This suggested the use of the approximation in (28c) (— s • e* ps |s • e es |). 

We wanted to check if the error indicators E*'s could provide a rough idea about the 
error level of the recovered solution. We also wanted to check if the local evaluations, 
at element level, of the E* could roughly describe the distribution of the error of the 
recovered solution. As far as we understand, there is no mathematical relation between 
E* and the error indicators E* that could make them provide accurate evaluations of 
E*. However, the numerical results obtained show that the E*'s capture the order of 
magnitude of E*. In particular, E% provides an accurate evaluation of E* ( E% ps E*) and 
very similar error distributions. 



6 Numerical Examples 

Quadrilateral elements have been considered in the examples presented in this section. 
6.1 Westergaard problem. XFEM solution 

Our first results were obtained as part of the development of an error estimator for XFEM. 
Let us consider the Westergaard problem corresponding to an infinite plate loaded at 
infinity with biaxial tractions <j xoo = cr yoo = and shear traction r^, presenting a 
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crack of length 2a as shown in Figure 1. The problem is solved using an enriched finite 
element approximation, combining the externally applied loads to obtain different loading 
conditions: pure mode I, pure mode II or mixed mode. 



Mi i a i 1 1 T 




1 J \ 1 \ a! \ I 



Figure 1: Example 1. Westergaard problem. Infinite plate with a crack of length 2a under 
uniform tractions (biaxial) and t^. Finite portion of the domain Qq, modelled with 
FE. 



The numerical model corresponds to a finite portion of the domain (a = 1 and b = 4 in 
Figure 1). The applied projected stresses for mode I are evaluated from the analytical 
Westergaard solution [6]: 



(To 



x cos v sin — + v — o m sin n cos 

2 y 2 J y \t\ 2 V 2 2 



x cos v sin — — v — o m sin n cos 

2 y 2 / y |t| 2 V 2 2 



2/ — 5 7F=F I m cos 7T + n sin ^ 

2 2 



(30) 
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and for mode II: 



2 ( y cos — + x sin — 



J/7712 



m cos — hn sin — 

tr \ 2 2 



a 2 r n 



y — s — = I m cos — \- n sin 
2 2 



(31) 



x cos y sin — + y — ^ man n cos 

2 2 J y \t\ 2 V 2 2 



where the stress fields are expressed as a function of x and y, with origin at the centre of 
the crack. The parameters t, m, n and are defined as 

t = (x + iy) 2 — a 2 = (x 2 — y 2 — a 2 ) + i(2xy) = m + in 

m = Rett) = Re(z 2 - a 2 ) = x 2 - y 2 - a 2 

W V ; (32) 

n = Im(t) = (z — a ) = 2xy 

<p = Arg(t) = Arg(m — in) with <\> E [— ir, it) , i 2 = — 1 



For the problem analysed, the exact value of the SIF is given by 

Ki, ex = cr^y/Tra K H ^ ex = r^^/iia (33) 



Material parameters are Young's modulus E = 10 and Poisson's ratio v = 0.333. We 
consider loading conditions in pure mode I with = 100 and = 0, pure mode II with 
fjoo = and Too = 100, and mixed mode with = 100 and r M = 100. In the numerical 
analyses, we use a geometrical enrichment defined by a circular fixed enrichment area 
B(xq, r e ) with radius r e = 0.5, with its centre at the crack tip Xq as proposed in [7]. For the 
extraction of the SIF we define a plateau function with radius r q = 0.9. Bilinear elements 
are considered in the models. For the numerical integration of standard elements we use 
a 2 x 2 Gaussian quadrature rule. The elements intersected by the crack are split into 
triangular integration subdomains that do not contain the crack. We use 7 Gauss points 
in each triangular subdomain, and a 5 x 5 quasipolar integration in the subdomains of the 
element containing the crack tip. We do not consider correction for blending elements. 
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Note that for this problem solved using XFEM, the displacements recovery technique 
was not implemented, therefore we used a stress recovery technique called SPR-C, see [ ] 
(this technique enforces the recovered stress field to satisfy at each patch the internal and 
boundary equilibrium equations and the compatibility equation. The recovery procedure 
also splits the stress field into singular and smooth parts: cr = c smo + cr sing , and then uses 
different recovery procedures for each part as described in [8] and [3]. Therefore, we used 
the estimation of the error in the displacements field given in (13), and used the following 
expression (equivalent to equation (27)) to evaluate an upper bound of the error of the 
recovered solution for a sequence of N meshes 



"(Oct I 



< e 



< 



(0*1 L 2 l h WlL 2 ~ l e WlL 2 H0Il 2 



U 



(AT) 



U 



«Il 2 I s «Il 2 



i (. i ) es \L 2 | S (i) 



L 2 



E, 



{i)UB 



N-l 



(34) 



then, we evaluated |e(7v) e s| L by extrapolation from |e(j_x). 




The error in energy norm of the recovered solution was evaluated using the following 
expression (similar to (28c)) where the term associated to the lack of boundary equilibrium 
was neglected because rs;0: 



E (i)3= / | s ' ( n (N) -u (i))| dQ = / \ s (i)- e (i)es\dtt; i = l,...,N-l, (35) 
Jn Jn 

we evaluated E* N ^ 3 by extrapolation from E? N _^ 3 . 

Figure 2 shows the evolution of the global effectivity index for the error estimate E 3 and 
the upper bound E^ B . We can observe that although E^ B bounds the exact error, it does 
not asymptotically converge to 1. We believe that we obtain bounds of the exact error 
because of the use of the CS inequality. Nevertheless, we can observe that the effectivity 
index of the error estimator in energy norm for the recovered solution E% captures the 
order of magnitude of the exact error, providing effectivities very close to unity. In some 
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cases, as in the top left graph, we can obtain worse accuracies in the last mesh because 
of the extrapolation procedure used. 
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Figure 2: Example 1. Westergaard problem. Effectivity for the estimate of the error of 
the recovered solution in mode I, mode II and mix mode. 



6.2 Square plate with 4 order displacements 

The FE software used with this and the following problems is based on the use of Cartesian 
grids independent of the shape of the problem to be analyzed. Elements fully located in 
the interior of the domain are considered as standard elements. Elements fully outside 
of the domain are not assembled. Finally, in elements cut by the boundary, element 
matrices are evaluated considering only the portion of the elements within the domain. 
The Lagrange multipliers technique has been used to impose a weak satisfaction of the 
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displacement constraints. The convergence rate of the error obtained is the theoretical rate 
expected for the standard FEM. The FE code also uses a hierarchical data structure to 
perform /i-adaptive analyses based on element splitting. Multi-point constraints are used 
to enforce C° continuity of the solution between adjacent elements of different refinement 
levels. This problem considers a A th order displacements field over an infinite domain. 
A 2x2 square portion has been modelled. The corresponding body loads and Neumann 
conditions have been imposed. A model of the problem is represented in Figure 3. 

u x (x, y) = x 4 + 5(x 3 y) — 3(x 2 y 2 ) + x 3 
u y (x, y) =y A - 6(y 2 x 2 ) + 3(yx 3 ) + 2y 

h ( ~3£ 
x[X,V) ~ 2(l + u)(2u-ly■■ 
(9x 2 - 12xy + Ay 2 - Ax + ... 
v[Ax 2 + 20x?/ - Ay 2 + Ax)) 

b * ix > y) = 2(l + v)(2v-iy 
(Ay 2 - 3x 2 + 2xy + is(8x 2 - 12xy) 
E = 1000 v = 0.3 

Figure 3: Example 2: 2 x 2 square model and analytical solution 

In Figure 4 we represent the evolution of the global effect ivity index 6 considering the 
error estimates E* defined in (28). Note that in this case we have considered the lack of 
equilibrium both in the internal and boundary equilibrium equations. 

The results show that the error estimators El, E\ and E\ capture the order of magnitude 
of the exact error in energy norm of the recovered solution. In particular, the results 
obtained with E\ provide error effect ivity indexes very close to one, the desired value. 
Taking into account the procedure described to obtain these error estimators, the results 
obtained, especially with E\, are surprisingly accurate. 

Moreover, in Figure 5 we show the local error evaluated by E%\k and the exact error of the 

1 1 2 

recovered field ||e*|| fe for a sequence of /i-adapted meshes, k indicates that the quantities 
are evaluated at element level. For this problem both results are quite similar, thus E^\k 
is a good indicator at local level of the error of the recovered solution. 
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e e 




dof dof 
Figure 4: Example 1. 2 x 2 square. Effect ivity index. Q4 (left) and Q8 (right). 

6.3 Pipe under internal pressure 

Let us now consider a pipe under internal pressure, as described in Figure 6. Figure 7 
shows the evolution of the effectivity index for the different error estimates. 

For this problem we have also evaluated the distribution of the exact energy norm of 
the recovered solution ||e*||^ and the distribution of the error estimate E£\k (subindex 
k indicates that these errors are evaluated locally, at elements). In Figure 8 we have 
represented these results for a sequence of /i-adapted meshes. We can observe that both 
error distributions are quite similar. These results show that the error estimator for the 
recovered solutions E$ has a good behavior at global level, but also at local level. 

6.4 L-Shape. Singular problem 

Consider the problem of an infinite plate with a V-notch subjected to tractions. We have 
considered the Mode I loading condition. The model of the problem, as defined in Figure 
9 has a singular solution in stresses. In this case we get similar results for the evolution of 
the global effectivity index for estimators E% and E£ (see Figure 10), although E% exhibits 
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Figure 5: Example 2. Square. Q4. Local exact error of the recovered solution (left), local 
error estimates using E% estimator (right). 
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Figure 6: Example 3: Pipe model and analytical solution 
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Figure 7: Example 3. Pipe. Effect ivity index. Q4 (left) and Q8 (right). 
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2G 
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r L [(k -Q{L+ 1)) cos (L$) - Lcos(((L - 2) cj>)] 



r L [(k + Q{L + 1)) sin (L4>) + Lcos{((L - 2) 4>)] 

2Gr 

£ = 1000 v = 0.3 

L = 0.544483736782464 Q = 0.543075578836737 
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Figure 9: Example 4: L-Shape model and analytical solution 
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Figure 10: Example 4: L-Shape. Effectivity index. Q4 (left) and Q8 (right). 



the best results for linear elements. In the last meshes of the analysis we always obtain an 
increase of the effectivity index. This is because in the first 3 meshes of each analysis the 
mesh increases its density around the reentrant corner increasing the refinement level as 
we get closer to the singular point. The FE code can reach up to 20 refinement levels. In 
the last mesh the refinement level around the singularity would require higher refinement 
levels. The result is that, as higher refinement levels cannot be reached, we obtain an 
area around the singularity with elements of uniform size as opposite to graded meshes 
towards the singularity. This produces pollution errors and a decrease in the accuracy of 
the error estimation that leads to worse effectivity indexes. 
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We show the local error index E^\ k and the exact error of the recovered solution, ||e*|| fc , 
in Figure 11. Observe that even for this singular problem the results are quite similar, 
concluding that, in this case E^\ k is a good error indicator. 

6.5 h-adaptive process 

In previous sections we have shown that the method presented for the estimation of error 
of the recovered solution, using the estimator E£, provides very accurate results with a 
very good global effectivity index. Now, we are going to show that this error estimation 
can be used in /i-adaptive processes.. 

During a /i-adaptive refinement process, using a standard FE compilation, the process 
reads as follows: 

1. Generate a FE mesh. 

2. Solve the FE problem. 

3. Estimate the error of the FE solution (locally and globally). 

4. If target error is smaller than the estimated error of the FE solution continue to 
step 5 else stop the process. 

5. Generate an /i-adapted FE mesh using the local FE error estimation 

6. Go to step 2. 

In this classical situation we are estimating the error of the raw FE solution, then we 
are using the raw FE solution (u h ,cr h ) as output. However, when we use our recovery 
procedure we have already available an improved solution (u*, <x*). So far, we were unable 
to estimate the error of this last solution, this output was not reliable. However, with 
the contribution presented in this report we have a way to obtain an accurate estimation 
of the error of this recovered solution: E%, thus we can use (u*, cr*) as the output of the 
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analysis. The information about the error estimation in the recovered solution could be 
used in the /i-adaptive refinement process to obtain a solution (u*, cr*) with the required 
accuracy: 

1. Generate a FE mesh. 

2. Solve the FE problem. 

3. Evaluate the local error estimate of FE solution (u h , <r h ). 

4. Evaluate the global value of the error estimate of the recovered solution (u*,cr*) 
(costless procedure). 

5. If target error is smaller than the error of the recovered solution continue to step 5 
else stop the process. 

6. Generate an /i-adapted FE mesh using the local FE error estimation 

7. Go to step 2. 

The /i-adaptive procedure is guided in both cases by the well established techniques based 
on the error estimation of the FE solution, which will simultaneously decrease the error 
of the FE solution and the error of the recovered solution. Note that the main difference 
between both processes is, simply, the stopping criterion. In this new approach the process 
stops when the estimated error of the recovered solution is smalled than the target error. 
The error in energy norm of the recovered solution is, in general, lower than the error in the 
FE solution and the convergence rate of the error in energy norm of the recovered solution 
is higher than that of the FE solution, thus, the recovered solution reaches the prescribed 
error level with less degrees of freedom than the FE solution. Therefore, this second 
method produces important savings in total computational cost of the analysis. Figures 
12 and 13 show the evolution of the errors during the /i-adaptive refinement process and 
the computational cost needed to obtain a certain accuracy level, for Example 3 (cylinder 
problem) for Q4 and Q8 elements. 

The horizontal black lines in figures 12 and 13 represent the error level of the solution 
prescribed by the analyst. Red and brown lines represent the error (exact and estimated) 
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(a) Convergence analysis. Relative error in (b) Computational cost analysis. Relative 
energy norm versus degree of freedom error in energy norm versus time. 

Figure 12: Example 3. /i-adaptive analysis with Q4 elements. The black line represents 
the prescribed relative error in energy norm (1%). 
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Figure 13: Example 3. /i-adaptive analysis with Q8 elements. The black line represents 
the prescribed relative error in energy norm (0.05%). 
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of the standard FE output (u h ,cr h ). Blue and green lines represent the error (exact and 
estimated) of the recovered field (u*,er*). We can observe that the error estimation of 
recovered solution accurately represents the exact error for both linear (Q4) and quadratic 
(Q8) elements. 

The figures show a considerable improvement of the /i-adaptive process to reach a pre- 
scribed error level. Note that, for Q4 elements, to reach the prescribed relative error in 
energy norm (1%), the standard /i-adaptive strategy based on the accuracy of the FE 
solution would stop the process after an analysis with 30492 dofs and 170.2 s whereas 
with the proposed /i-adaptive method based on the accuracy of the recovered solution the 
process would stop with only 654 dofs and 11.7 s. With Q8 elements reaching a 0.05% 
prescribed error, the standard /i-adaptive procedure would stop the process after an anal- 
ysis wiht 41724 dofs and 54.4 s, whereas the proposed method would only require 3728 
dofs and 21.8 s. 



7 Conclusions 



Error estimation of the recovered solution could be very useful to decrease the compu- 
tational cost of FE analysis. Error estimator efficiently predicts the error in energy 
norm of the recovered solution both at local and global levels. This accuracy is not math- 
ematically justified. The results obtained definitively require the support of mathematical 
proofs. 
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